from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from numpy import arange

plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

param = np.array([
    [30, 50, 6, 18.193],
    [30, 70, 8, 9.12],
    [30, 90, 10, 9.276],
    [40, 50, 6, 30.243],
    [40, 70, 8, 23.856],
    [40, 90, 10, 35.673],
    [50, 50, 6, 6.586],
    [50, 70, 8, 6.586],
    [50, 90, 10, 8.323],
])


def f(T, H, SC, k_1, k_3, k_4, k_5, k_6, k_H, S_C, S_S, eta0, n, alpha):
    # 变量
    solid_content = SC / 100  # 固含量
    T = T + 273.15  # 温度，单位为开尔文

    # 常量
    C_0 = 24e-3  # 比例系数，用于确定DMF质量与时间的关系
    m_D_0 = 24e-3  # 初始时刻DMF的质量，单位为kg
    m_S = 6e-3  # 环丁砜总的质量，单位为kg。无论是否溶解，环丁砜不会从体系中消失，总的质量是不变的
    m_C = (m_D_0 + m_S) * solid_content  # 初始时刻醋酸纤维素的质量，单位为kg，与初始固含量有关
    ro_D = 0.948e3  # DMF的密度，单位为kg/m^3
    ro_S = 1.261e3  # 环丁砜的密度，单位为kg/m^3
    ro_C = 1.3e3  # 醋酸纤维素的密度，单位为kg/m^3
    V = m_D_0 / ro_D + m_S / ro_S + m_C / ro_C  # 总体积，单位为m^3，认为不变
    A = 6.09451
    B = 2725.96
    C = 28.209  # 由温度计算饱和汽压时的三个常数
    A_0 = -k_1 * (10 ** (A - B / (T + C))) * 133.322 * (1 - k_H * H / 100) / ro_D / V
    n_correction_factor = n  # 斯托克斯 - 爱因斯坦方程中的修正系数，对于大分子乃至宏观颗粒是6，对于小分子则体积越小，该值也会越小
    k = 1.380649e-23  # 玻尔兹曼常数，单位为J/K
    r = 0.3413e-9  # 环丁砜的分子半径
    N_A = 6.02214076e23  # 阿伏伽德罗常量，单位为mol^-1

    # 函数
    # DMF在蒸发的过程中，其质量随时间变化的函数
    def m_D(t):
        return C_0 * np.exp(A_0 * t)

    # 未能溶解而析出成液滴的环丁砜的质量随时间变化的函数
    def m_S_out(t):
        return max(0, m_S - m_D(t) * S_S)

    # 未能溶解而析出成固体的醋酸纤维素的质量随时间变化的函数
    def m_C_out(t):
        return max(0, m_C - m_D(t) * S_C)

    # 析出的醋酸纤维素的体积占比随时间变化的函数
    def phi_C_out(t):
        return m_C_out(t) / ro_C / V

    # 不同温度下混合溶液总的粘度随时间变化的函数
    def eta(temp, t):
        return eta0 * (1 + 2.5 * phi_C_out(t))

    # 扩散系数。假定一个小液滴内仅有一个分子，且分子为球形
    def D(temp, t):
        return k * temp / (n_correction_factor * np.pi * r * eta(temp, t))

    # 液滴运动的平均速度
    def v(temp, t):
        return k_3 * (D(temp, t) ** 0.5)

    # 析出的环丁砜的分子数密度随时间变化的函数
    def n_molecular_number_density(t):
        return m_S_out(t) / (ro_S * 4 / 3 * np.pi * r ** 3) / V

    # 小液滴间的碰撞频率
    def Z(temp, t):
        return (2 ** 0.5) * n_molecular_number_density(t) * np.pi * (r ** 2) * v(temp, t)

    tf = np.log(alpha) / A_0  # 对于每种情况，在反应了一定时间后停止

    dy = lambda m_s, t: k_4 * Z(T, t) + k_5 * v(T, t) * (m_S_out(t) - m_s) / V * m_s
    t = arange(1, tf, 0.01)
    sol = odeint(dy, 0, t)

    if len(sol.T[0]) == 0:
        re = 10000
    else:
        re = k_6 * sol.T[0][-1]
    return re


def pred(k_1, k_3, k_4, k_5, k_6, k_H, S_C, S_S, eta0, n, alpha):
    predictions = []
    for i in range(param.shape[0]):
        T = param[i][0]
        H = param[i][1]
        SC = param[i][2]
        re = f(T, H, SC, k_1, k_3, k_4, k_5, k_6, k_H, S_C, S_S, eta0, n, alpha)
        predictions.append(re)
    return predictions


def target():
    targets = []
    for i in range(param.shape[0]):
        t = param[i][-1]
        targets.append(t)
    return targets


def loss(pred, target):
    loss = 0
    for i in range(param.shape[0]):
        loss += np.linalg.norm(pred[i] - target[i])
    return loss


# 真实值和预测值
y_true = np.array([])
y_pred = np.array([])


def printre(kkb):
    k_1, k_3, k_4, k_5, k_6, k_H, S_C, S_S, eta0, n, alpha = kkb
    t = target()
    global y_true, y_pred

    for i in range(param.shape[0]):
        T = param[i][0]
        H = param[i][1]
        SC = param[i][2]
        re = f(T, H, SC, k_1, k_3, k_4, k_5, k_6, k_H, S_C, S_S, eta0, n, alpha)
        tt = t[i]
        print(f"预测值{re},真实值{tt}")
        y_true = np.append(y_true, [tt])
        y_pred = np.append(y_pred, re)


def plot_SA_history():
    plt.figure(figsize=(14, 10))

    # 子图1：损失函数变化
    plt.subplot(2, 2, 1)
    plt.semilogy(history['best_loss'], 'r-', label='Best Loss')
    plt.semilogy(history['current_loss'], 'b--', alpha=0.5, label='Current Loss')
    plt.xlabel('Iteration')
    plt.ylabel('Loss')
    plt.title('Loss变化曲线')
    plt.legend()
    plt.grid(True, which="both", ls="--")

    # 子图2：温度衰减曲线
    plt.subplot(2, 2, 2)
    plt.plot(history['temperature'], 'g-')
    plt.xlabel('Iteration')
    plt.ylabel('Temperature')
    plt.title('温度下降曲线')
    plt.grid(True, ls="--")

    # 子图3：k参数演化
    plt.subplot(2, 2, 3)
    k_array = np.array(history['k_params'])
    for i in range(5):
        plt.plot(k_array[:, i], label=f'k_{i + 1}')
    plt.xlabel('Iteration')
    plt.ylabel('k参数值')
    plt.title('k参数变化曲线')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, ls="--")

    # 子图4：kh参数演化
    plt.subplot(2, 2, 4)
    kh_array = np.array(history['kh_params'])
    label_list = ["醋酸纤维素溶解度", "环丁砜溶解度", "斯托克斯-爱因斯坦方程修正系数", "反应比例"]
    for i in range(4):
        plt.plot(kh_array[:, i], label=label_list[i])
    plt.xlabel('Iteration')
    plt.ylabel('物理参数值')
    plt.title('物理参数变化曲线')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, ls="--")

    plt.tight_layout()
    plt.show()


# 在SA函数外部初始化记录容器
history = {
    'temperature': [],
    'best_loss': [],
    'current_loss': [],
    'k_params': [],
    'kh_params': []
}


def SA(t0, tf, alpha, iter):
    global history
    t = t0
    # k_10, k_30, k_40, k_50, k_60, k_H0, S_C0, S_S0, eta00, n0, param_alpha0 = np.array([1e-2, 1e+0, 1e-2, 1e-1, 1e-1, 1e-1, 1, 0.8, 0.8, 1, 0.2])
    k_10, k_30, k_40, k_50, k_60, k_H0, S_C0, S_S0, eta00, n0, param_alpha0 = np.array([0.002, 0.2, 0.2, 0.04, 0.2, 0.4, 0.7, 0.5, 0.7, 2.4, 0.3])
    k_1c, k_3c, k_4c, k_5c, k_6c, k_Hc, S_Cc, S_Sc, eta0c, nc, param_alphac = k_10, k_30, k_40, k_50, k_60, k_H0, S_C0, S_S0, eta00, n0, param_alpha0

    lc = loss(pred(k_1c, k_3c, k_4c, k_5c, k_6c, k_Hc, S_Cc, S_Sc, eta0c, nc, param_alphac), target())

    k_1b, k_3b, k_4b, k_5b, k_6b, k_Hb, S_Cb, S_Sb, eta0b, nb, param_alphab = k_1c, k_3c, k_4c, k_5c, k_6c, k_Hc, S_Cc, S_Sc, eta0c, nc, param_alphac
    lb = loss(pred(k_1b, k_3b, k_4b, k_5b, k_6b, k_Hb, S_Cb, S_Sb, eta0b, nb, param_alphab), target())
    for i in range(iter):
        k_1n = k_1c + np.random.normal(0, 1e-4)
        k_1n = np.clip(k_1n, 1e-4, 5e-3)
        k_3n = k_3c + np.random.normal(0, 1e-1)
        k_3n = np.clip(k_3n, 1e-1, 1e+0)
        k_4n = k_4c + np.random.normal(0, 1e-1)
        k_4n = np.clip(k_4n, 1e-1, 1e-0)
        k_5n = k_5c + np.random.normal(0, 1e-3)
        k_5n = np.clip(k_5n, 1e-2, 1e-1)
        k_6n = k_6c + np.random.normal(0, 1e-2)
        k_6n = np.clip(k_6n, 1e-1, 1e+0)
        k_Hn = k_Hc + np.random.normal(0, 1e-2)
        k_Hn = np.clip(k_Hn, 1e-1, 1e-0)
        S_Cn = S_Cc + np.random.normal(0, 1e-1)
        S_Cn = np.clip(S_Cn, 1e-1, 1e+1)
        S_Sn = S_Sc + np.random.normal(0, 0.1)
        S_Sn = np.clip(S_Sn, 0, 1)
        eta0n = eta0c + np.random.normal(0, 0.1)
        eta0n = np.clip(eta0n, 2, 3)
        nn = nc + np.random.normal(0, 0.1)
        nn = np.clip(nn, 0, 6)
        param_alphan = param_alphac + np.random.normal(0, 0.01)
        param_alphan = np.clip(param_alphan, 0.1, 1)
        '''
        [0.014840809916094374, 0.1430051145124099, 0.07345543679797854, 0.10474345856177354, 0.35796509420390066, 0.39955637570468105, 0.835754321687266, 0.8707603453814515, 0.7408325702069575, 5.551160191461859, 0.001]
        预测值17.975241480357127,真实值18.193
        预测值19.354695773656555,真实值9.12
        预测值20.968373550688504,真实值9.276
        预测值11.136179654426675,真实值30.243
        预测值12.118839124052068,真实值23.856
        预测值13.3216971785971,真实值35.673
        预测值6.834646178989657,真实值6.586
        预测值7.477928380971328,真实值6.586
        预测值8.259533902404392,真实值8.323
        MAPE: 58.15%
        
        [0.0010660234370888787, 0.34052550181423463, 0.11440264521150946, 0.03446581365646933, 0.30328518806514704, 0.32248515839275166, 0.9118799744432465, 0.6509464523758399, 0.700032233700628, 2.495457103926549, 0.2574903993186905]
        预测值17.973155027875798,真实值18.193
        预测值19.167187958646682,真实值9.12
        预测值20.507658638913224,真实值9.276
        预测值11.139508097625269,真实值30.243
        预测值11.999026665286273,真实值23.856
        预测值13.003991335670438,真实值35.673
        预测值6.853522977744103,真实值6.586
        预测值7.402026716535338,真实值6.586
        预测值8.041227613622853,真实值8.323
        MAPE: 47.63%
        
        [0.0001, 0.04062609086637083, 0.20576415942501514, 0.1656558178892305, 0.2777948424758192, 0.33712968164168544, 0.7454938345585671, 0.8711713046463835, 0.8259010553802277, 0.70865543089256, 0.2308007539262029]
        预测值13.939596208660502,真实值18.193
        预测值14.307984574702267,真实值9.12
        预测值14.657446160466122,真实值9.276
        预测值10.719022184775255,真实值30.243
        预测值11.29540231580299,真实值23.856
        预测值11.902760260729094,真实值35.673
        预测值7.324940698412884,真实值6.586
        预测值7.857626388089587,真实值6.586
        预测值8.457383068512778,真实值8.323
        MAPE: 39.36%
        
        [0.010789226098564548, 1.387920573590249, 0.019230368698797965, 0.0767177933026201, 0.2405455765172372, 0.36313831619148434, 0.7716201213393014, 0.7747745286247333, 0.8833722048753654, 0.783226713389143, 0.05154546949979516]
        预测值12.36122596728503,真实值18.193
        预测值12.53378203807114,真实值9.12
        预测值12.695530089500505,真实值9.276
        预测值10.291385127859856,真实值30.243
        预测值10.694381580546208,真实值23.856
        预测值11.114930136842254,真实值35.673
        预测值7.445403721796783,真实值6.586
        预测值7.905239770173509,真实值6.586
        预测值8.432947916742991,真实值8.323
        MAPE: 36.75%
        
        [0.004897901963678431, 1.0838414870419346, 0.021763427795672512, 0.10706129271564353, 0.2846622679320061, 0.33190409579994723, 0.8972423866200637, 0.9187726268393154, 0.9266555816222484, 1.038394396980951, 0.07974238424869429]
        预测值12.069769343095455,真实值18.193
        预测值12.21850632162352,真实值9.12
        预测值12.34422324062141,真实值9.276
        预测值10.169133962092493,真实值30.243
        预测值10.54194234615023,真实值23.856
        预测值10.882726344065869,真实值35.673
        预测值7.439326242415858,真实值6.586
        预测值7.86286990044933,真实值6.586
        预测值8.321343186112763,真实值8.323
        MAPE: 36.08%
        '''

        pre = pred(k_1n, k_3n, k_4n, k_5n, k_6n, k_Hn, S_Cn, S_Sn, eta0n, nn, param_alphan)
        ln = loss(pre, target())

        if ln < lc or np.random.rand() < np.exp(-(ln - lc) / t):
            k_1c, k_3c, k_4c, k_5c, k_6c, k_Hc, S_Cc, S_Sc, eta0c, nc, param_alphac = k_1n, k_3n, k_4n, k_5n, k_6n, k_Hn, S_Cn, S_Sn, eta0n, nn, param_alphan
            lc = ln
            if lc < lb:
                k_1b, k_3b, k_4b, k_5b, k_6b, k_Hb, S_Cb, S_Sb, eta0b, nb, param_alphab = k_1c, k_3c, k_4c, k_5c, k_6c, k_Hc, S_Cc, S_Sc, eta0c, nc, param_alphac
                lb = lc
        kkb = [k_1b, k_3b, k_4b, k_5b, k_6b, k_Hb, S_Cb, S_Sb, eta0b, nb, param_alphab]
        t *= alpha
        if t < tf:
            break
        kc = np.array([k_1c, k_3c, k_4c, k_5c, k_6c, k_Hc, S_Cc, S_Sc, eta0c, nc, param_alphac])
        print(f"iter{i}, lb{lb},kc{kc},ln{ln},")
        k_params = np.array([k_1c, k_3c, k_4c, k_5c, k_6c])
        kh_params = np.array([k_Hc, S_Cc, S_Sc, eta0c])
        history['temperature'].append(t)
        history['best_loss'].append(lb)
        history['current_loss'].append(lc)
        history['k_params'].append(k_params.copy())
        history['kh_params'].append(kh_params.copy())
    return kkb


re = SA(300, 0.001, 0.9, 1000)
print(re)
printre(re)
plot_SA_history()

# 计算 MAPE
mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
print(f"MAPE: {mape:.2f}%")
